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ABSTRACT 


The  work  described  is  a  continuation  of  that  reported 
previously  in  AEDC-TN-6I-65,  AEDO-TDR-62-131,  and  AEDC-TDR-63-82. 
The  present  report  contains  a  general  method  for  Integrating 
numericali.y  through  a  saddle-point  singularity  of  an  ordinary 
differential  equation.  The  solution  of  the  differential  equation 
is  assumed  to  depend  also  on  the  value  of  a  parameter,  such  as 
the  mass  flow.  The  method  is  thus  applicable  to  a  wide  assort¬ 
ment  of  gas-dynamic  problems  including  one-dimensional  nonequi¬ 
librium  nozzle  flow  and  two-phase  nozzle  flow. 

Specific  application  is  made  to  nonequilibrium  nozzle  flow, 
and  the  results  of  this  application  are  presented  and  discussed. 
The  method  proved  to  be  numerically  accurate  without  requiring  an 
exceedingly  precise  estimate  for  the  critical  mass  flow. 

The  work  also  includes  a  modification  of  the  method  for 
calculating  approximate  equilibrium  nozzle  flows  first  given  i n 
AEDC-TDR-62-131. 
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1.0  INTRODUCTION 

This  Is  the  fourth  and  final  report  In  a  series  concerning  the 
calculation  of  the  reacting  flow  of  a  complex  gas  in  the  nozzle  of  a 
hypersonic  wind  tunnel.  The  preceding  three  reports,  hereafter 
referred  to  as  Parts  I,  II,  and  III,  dealt  with  the  following  topics: 

(a)  Part  I,  hy  Vincent!  (I961),  descrihes  a  five-species  model 
for  air,  governed  hy  eight  chemical  kinetic  reactions.  A 
method  is  given  for  the  numerical  calculation  of  the  one- 
dlmenslonal  nonequilibrium  flow  of  this  gas  through  a 
hypersonic  nozzle.  A  specific  calculation  carried  out  on 
a  high-speed  digital  computer  revealed  that  the  method 
required  too  much  computer  time  to  be  practical  for 
engineering  purposes. 

(b)  Part  II,  by  Emanuel  and  Vincent!  (196?^  dcsei’-! -  'ethod 

...  rHinibrlntf  ^  pie, 

...oxve  amount  of 

machine  time.  This  1  poio  u,j.su  contains  an  approximate 
but  X’elatlvely  simple  method  for  calculating  equilibrium 
nozzle  flows. 

(c)  Part  III,  by  Emanuel  (1963),  contains  the  analytical  basis 
for  the  numerical  method  described  In  Part  II.  The  report 
Is  thus  concerned  primarily  with  the  interaction  of  a 
"stiff",  ordinary  differential  equation  and  various  inte¬ 
gration  procedures. 

The  method  described  in  Part  II  requires  that  the  forward  integra¬ 
tion  of  the  nonequilibrium  equations  proceed  from  an  initial  point  that 
is  in  equilibrium.  This  point  is  chosen  somewhat  upstream  of  the 
nozzle  location  at  which  the  chemistry  first  departs  appreciably  from 
the  equilibrium-flow  solution.  For  a  sufficiently  high  stagnation 
pressure,  this  point  may  be  taken  in  the  supersonic  portion  of  the  nozzle. 
In  this  situation  the  nonequilibrium  flow  is  readily  calculated  by  the 
method  given  in  Part  II.  For  low  stagnation  pressures,  however,  the 
chemistry  first  departs  from  the  equilibrium-flow  solution  at  a  subsonic 
nozzle  location,  and  the  forward  integration  of  the  equations  must  now 
proceed  from  a  subsonic  initial  point.  In  this  case  the  method  given 
in  Part  II  must  be  modified  to  incorporate  a  procedure  for  dealing  with 
the  singularity  that  occurs  at  the  sonic  point.  This  report  describes 
in  detail  one  such  method. 

The  sonic-point  singularity  occurs  in  both  frozen  and  nonequilibrium 
nozzle  flows.  In  both  instances  it  is  a  saddle  point  that  occurs  when 
the  numerator  and  denominator  of  the  right-hand  side  of  an  ordinary  dif¬ 
ferential  equation  simultaneously  become  zero.  The  differential  equation 
involved  is  a  gas -dynamic  one.  In  the  frozen-flow  case,  for  example, 
it  is  given  by 
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dx 


(1-1) 


where  T  Is  the  temperature,  x  the  distance  along  the  nozzle  axis, 
y  the  ratio  of  specific  heats,  Mj.  the  Mach  number  based  on  the  frozen 
speed  of  sound,  and  A  the  nozzle  area.  For  frozen  flow,  the  location 
of  the  singularity  is  at  the  throat  and  the  critical  mass  flow  is  easily 
determined.  Eijuatlon  (l-l)  assumes  that  the  nozzle  area  is  a  given 
function  of  x  ,  as  is  generally  the  ease,  while  the  temperature  is 
the  unknown  dependent  variable.  If  the  role  of  the  two  variables  is 
interchanged,  J.e.,  the  temperature  is  a  given  function  of  x  and  the 
area  is  the  unknown  dependent  variable,  then  the  singularity  is 
removed.  This  result  is  evident  when  equation  (l-l)  is  rewritten  as 
follows ; 


Thus,  the  singularity  in  one -dimensional  nozzle  flow  occurs  only  when 
the  doubled-valued  variable  A  is  taken  as  the  given  function  of  x  . 

In  nonequilibrium  nozzle  flow  the  situation  is  more  complicated. 

In  particular,  the  location  of  the  saddle  point  and  the  critical  mass 
flow  are  not  known  a  priori.  It  is  therefore  necessary  to  guess  a 
value  for  the  mass  flow  before  starting  the  numerical  integration. 

This  value,  of  course,  will  probably  differ  from  the  critical  one,  and 
Lhe  resulting  solution  will  either  be  subsonic  for  the  entire  nozzle, 
or  an  infinity,  which  terminates  the  solution,  will  occur  upstream  of 
the  singularity. 

The  problem  dealt  with  in  this  report  may  be  stated  as  follows: 

A  solution  is  to  be  obtained  that  passes  smoothly  through  the  saddle- 
point  singularity  of  an  ordinary  differential  equation.  This  solution 
furthermore  depends  on  a  parameter,  such  as  the  mass  flow,  whose  precise 
value  is  unknown.  As  such,  this  problem  is  encountered  not  only  in 
computing  nonequilibrium  nozzle  or  diffuser  flow  but  occurs  also  in  two- 
phase  nozzle  or  diffuser  flow  (see  Glauz  (1962)).  It  also  appears  in 
Inviscld  flow  about  a  two-dimensional  or  axisymmetric  blunt  body  in  a 
supersonic  stream  when  such  flow  is  computed  by  the  method  of  Integral 
relations  (see  Belotserkovskii  (l95T)  and  (1958)  on  Xerikos  and 
Anderson  (1962)). 

A  general  solution  to  the  foregoing  problem  is  outlined  in 
Chapter  2,  This  solution  takes  the  form  of  a  well-defined  numerical 
method,  hereafter  referred  to  as  the  optimum-point  method,  that  is  well 
suited  for  use  with  high-speed  digital  computers.  The  method  has  the 
important  advantages  of  being  relatively  simple  to  apply,  not  requiring 
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an  exceedingly  precise  estimate  of  the  parameter,  and  yielding  a  numeri¬ 
cally  accurate  solution.  Specific  application,  described  in  Chapter  3 
and  in  Appendix  I,  was  made  to  nonequilibrium  nozzle  flows.  Numerous 
calculations  on  an  IBM  7090  computer  verified  the  simplicity  and  accuracy 
of  the  method  in  this  application.  The  results  of  one  of  these,  along 
with  the  corresponding  equilibrium  nozzle -flow  solution,  is  presented  in 
Chapter  U.  This  chapter  also  discusses  certain  gas-dynamic  results 
stemming  from  this  calculation.  Finally,  Appendix  II  modifies  the  approx¬ 
imate  equilibrium  method  first  presented  in  Part  II. 

The  singularity  problem  can  also  be  attacked  by  means  of  trlal-and- 
error  procedures,  wherein  numerous  Integrations  are  performed,  each  with 
a  slightly  different  value  for  the  mass  flow  (see,  for  example, 

Reichenbach  (I96O)).  Once  a  subsonic  solution  is  available  that  is 
sufficiently  close  to  the  critical  one,  the  integration  can  then  be 
started  from  a  supersonic  initial  point  obtained  by  extrapolating  the 
subsonic  solution  across  the  singularity.  This  approach,  however,  consumes 
a  considerable  amount  of  computer  time  and  generally  results  in  a  more 
accurate  estimate  for  the  critical  mass  flow  than  is  actually  needed. 

Other  methods  also  exist  for  dealing  with  saddle-point  singularities. 
These  methods  are  applicable  only  in  special  clrcumsbances.  For  example, 
Eschenroeder,  Boyer  and  Hall  (1962)  outline  a  method  for  dealing  with 
the  singularity  in  nonequilibrium  nozzle  flow.  Their  method  is  based  on 
the  assumption  that  the  nonequilibrium  value  for  the  density  differs  neg¬ 
ligibly  from  its  equilibrium-flow  value  in  the  subsonic  portion  of  the 
nozzle.  The  equilibrium  value  of  the  density  is  then  used  as  the  inde¬ 
pendent  variable  to  calculate  the  subsonic  portion  of  the  nozzle.  Another 
special  method  designed  for  two-phase  flows  is  due  to  Glauz  (1962).  This 
method  Introduces  a  new  dependent  variable  N  defined  by 

N  =  I  (M^  -1)^  .  (1-3) 

Transformation  (l-3)  does  remove  the  singularity  from  equation  (l-l). 

The  transformed  equation,  however,  must  now  simultaneously  satisfy  the 
interior-boundary  conditions  N  =  0  and  (dN/dx)  =  0  . 

The  author  wishes  to  express  his  sincere  gratitude  to  Professor 
Walter  G.  Vincentl  for  his  encouragement  and  assistance,  Thanks  are  also 
due  to  Mrs.  Lita  Emanuel  for  editorial  assistance  and  to  Messrs.  H.  L. 
Mitchell  and  Lyle  B.  Smith  for  computer  programming.  During  the  work 
the  author  was  supported  in  part  by  a  research  grant  to  Stanford  University 
from  the  National  Science  Foundation. 
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2.0  OPTIMUM-POINT  METHOD 

Before  the  optimum-point  method  is  described  in  its  simplest  form, 
some  pertinent  facts  concerning  saddle  points  will  be  presented. 
Consider  the  single  ordinary  differential  equation 


^  =  P(x,y) 

dx  Q{x,y) 


(2.1) 


with  the  Initial  condition  y(x  )  =  y  .  In  addition,  let  the  right- 
hand  side  of  equation  (2-l)  depend  on  a  parameter  a  not  explicitly 
shown  in  the  equation.  Consequently,  the  solution  of  equation  (2.1), 

y  =  y(x;a)  ,  (2-2) 


depends  on  o:  .  We  also  require  that  for  a  unique  value  of  a  ,  refer¬ 
red  to  as  the  critical  value  ,  the  differential  equation  has  a 

saddle  point  at  (x'’*',y*)  ,  as  shown  in  Plgure  1.  For  values  of  a 
different  from  ,  the  integral  curves  lie  either  above  or  below  the 

critical  solution  =  y(x;a  )  .  Only  two  solutions  of  equation  (2-1) 

actually  pass  through  the  saddle  point.  One  of  these  originates  at  a 
point  above  y*  ,  while  the  other  originates  below.  We  may  restrict  the 
discussion  to  the  case  yQ  >  y*  without  Toss  of  generality. 

In  general,  the  critical  value  of  a  is  not  known  a  priori.  In 
this  case,  for  some  value  of  0!  ,  say  (see  Figure  1),  the  integral 

curve  y(xja2^)  has  a  minimum  where  P(x,y(x;%))  is  zero.  This  zero 
generally  does  not  occur  at  x  =  x*  .  Similarly,  each  integral  curve 
below  the  critical  one  has  a  point  where  Q  is  zero  and  the  elope  is 
infinite.  This  zero  does  not  generally  occur  at  y  -  y*  .  Thus,  the 
exact  location  of  the  saddle  point  (x^jy*)  is  also  not  known  a  priori. 

It  is  well  known  that  the  critical  solution  in  a  neighborhood  of 
the  saddle  point  is  linear.  The  slope  of  the  critical  solution  at  the 
saddle  point  can  be  determined  if  L 'Hospital's  rule  can  be  applied  to 
the  indeterminate  form  P/Q  .  When  the  location  of  the  saddle  point  is 
not  known  a  priori,  however,  L'Hospital's  rule  cannot  be  applied 
directly. 

Finally,  we  note  that  if  (X  is  close  to  ,  then  the  Integral 

curve  y(x;a)  remains  close  to  y^^.  until  either  P  or  Q  is  nearly 
zero.  Furthermore,  the  slope  dy(x;a)/dx  also  remains  close  to 
dycr/dx  until  the  Integral  curve  y(x;a)  approaches  the  saddle  point. 
Thus,  if  a  is  sufficiently  close  to  ,  there  will  exist  a  point 

(x  pt,yqp^)  on  y  =  y(x;a)  ^  close  to  the  saddle  point,  where  the  slope 
dy(x;Q)7ax  best  approximates  the  slope  of  the  critical  solution  evaluated 
at  the  saddle  point.  The  function  of  the  optimum-point  method  is  to 
locate  this  point.  It  must  be  re-emphasized  that  the  location  of  the 
optimum- point  cannot  depend  explicitly  on  the  critical  value  ,  on 
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the  location  of  the  saddle  point  (x*,y*)  ,  or  on  the  slope  of  the  criti¬ 
cal  solution  at  the  saddle  point,  since  all  these  quantities  are  unknown. 

The  actual  location  of  the  optimum  point  depends  on  a  criterion  based 
on  L 'Hospital's  rule.  When  this  rule  is  applied  to  equation  (2-1),  we 
obtain 


dy  dx 
dx  “  ^  ’ 
dx 


(2-3) 


where  the  derivative  d(  )/dx 


is  given  by 


_d S  ^  dy  5 

dx  “  §x  dx  5y  * 


When  dy/dx  is  eliminated  from  equations  (2-1)  and  (2-3),  the  following 
result  is  obtained: 


--§-  =  0 
dx  dx 


(2-4) 


This  result  is  exact,  of  course,  only  at  the  saddle  point.  The  left-hand 
side  of  equation  (2-4),  which  has  the  same  dimensions  as  x  ,  in  general 
is  not  zero  when  evaluated  away  from  the  saddle  point.  A  new  parameter 
XpQ  is  therefore  defined  as 


■  -  Z. 

■PQ  "  ^  ^  ■ 

dx  dx 


(2-:^) 


This  parameter  is  evaluated  along  a  specific  Integral  curve  during  the 
forward  integration  of  equation  (2-l).  How  XpQ  is  evaluated  numerically 
is  discussed  in  the  next  chapter.  The  optimum  point  is  then  defined  as 
the  point  nearest  to  (x*,y*)  where  I^pqI  Is  e  minimum. 

We  now  indicate  why  one  of  the  several  minimum  values  of  I^^pqI 
that  may  occur  gives  the  point  at  which  the  slope  dy(xjO:)/dx  best 
approximates  the  slox)e  of  the  critical  solution  evaluated  at  the  saddle 
point.  We  first  assume  that  P  and  Q  are  continuously  differentiable, 
as  is  generally  the  case.  Next,  let  us  rev/rlte  Xpq  as  follows: 


X 


PQ 


dx 
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Since  P  and  Q  are  continuously  differentiable,  the  ratio 
(dP/dx)/(dQ/dx)  may  be  evaluated  at  the  saddle  point,  which  is  assumed 
to  be  close  to  the  optimum  point.  The  ijuantity  in  parentheses  is  then 
approximately  the  difference  between  the  slope  dy(x;Q;)/dx  and  the  slope 
of  the  critical  solution  at  the  saddle  point.  Therefore,  this  minimum 
value  of  I^pqI  corresponds  to  the  optimum  point. 

As  o:  tends  toward  ,  the  integral  curve  y  =  y(x;a)  tends 

toward  the  critical  solution.  Thus,  the  optimum  point  will  also  tend 
toward  the  saddle  point  as  shown  by  the  dots  in  Figure  1.  All  of  the  non- 
eq^uilibrium  nozzle -flow  calculations  carried  out  in  the  course  of  this 
work  are  in  accord  with  the  foregoing  conclusions. 

On  the  basis  of  the  foregoing  ideas,  an  approximate  numerical  solu¬ 
tion  that  passes  through  the  saddle  point  can  be  constructed  as  follows; 

(a)  The  numerical  integration  of  equation  (2-l),  for  a  particular 

value  of  a  ,  proceeds  from  until  the  optimum  point 

is  reached. 

(b)  For  values  of  x  in  an  interval  S  x  ^  x^^  ,  where  the 

transfer  point  x^^  is  defined  as  some  suitable  point  greater 
than  x*  ,  the  linear  approximation 


is  used. 

(c)  For  values  of  x  greater  than  x^^  equation  (2-1)  is  again 
integrated  numerically  with  the  initial  condition 

^Kr)  =  Kr-^opt)(S)^p^+  ^opt  •  (2-7) 

This  scheme  is  deficient  in  that  dy/dx  will  have  a  discontinuity  at 
x^j,  since  P/Q  ,  evaluated  at  (x^j.,y(x^j,) )  ,  will  not  generally  be  equal 
to  (dy/dx)Qp^  .  This  difficulty  can  be  avoided  if  it  is  possible  to 
rewrite  equation  (2-1)  in  an  alternative  form  that  does  not  contain  a 
singularity  in  the  region  S  x  ^  x,  .  In  this  case,  both  equations 

(2-1)  (in  its  alternative  form)  and  (2-6)  apply  and  consequently  p/q  , 
evaluated  at  (x^p,y(x^p))  ,  will  equal  (dy/dx)Qp^  and  no  discontinu¬ 
ity  occurs.  Chapter  3.0  demonstrates  how  this  is  accomplished  for 
nozzle-flow  calculations. 

The  foregoing  scheme  is  practical  only  if  |a-apj.|  is  small.  By 
small  we  mean  that  a  and  agree  to  three  or  four  signif icantf Igures. 

This  is,  in  fact,  a  considerable  improvement  over  other  saddle -point 
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methods,  which  frequently  require  much  greater  accuracy  for  the  value  of 
a  .  The  condition  that  (a-a^-p|  he  small  can  he  satisfied,  for  example, 
hy  requiring  that  both  |P(Xopt>yopt^  I  l*^(’'^opt'^opt)  I  small. 

The  numerical  calculation  can^then  he  fully  automated  hy  including  pro¬ 
vision  for  restarting  the  calculation  at  when  one  or  both  of  these 

conditions  is  not  satisfied.  The  new  calculation  would,  of  course,  pre¬ 
sumably  use  a  value  for  a  considerably  closer  to  than  the  preced¬ 

ing  one. 

An  additional  advantage  of  the  optimum-point  method  is  that  no  special 
choice  for  the  variables  is  necessary.  Furthermore,  those  equations  not 
containing  the  saddle-point  singularity  arc  not  affected. 


3.0  APPLICATION  TO  NONEQUILIBRIUM  NOZZLE  FLOW 
3.1  INTRODUCTORY  REMARKS 


This  chapter  applies  the  optimum-point  method  to  the  saddle-point 
singularity  that  occurs  in  steady  one-dimensional  nozzle  flow.  For  con¬ 
venience,  the  discussion  will  be  limited  to  nonoquilibrlum  nozzle  flow, 
although  the  method  is  equally  applicable  to  nonequilibrium  diffuser  flow 
and  two -phase  nozzle  or  diffuser  flow.  In  addition,  the  method  has  been 
applied  successfully  to  the  approximate  equilibrium  nozzle-flow  method 
given  in  Part  II. 

Specifically,  this  chapter  describes  how  the  nonequilibrium  nozzle - 
flow  equations  can  be  integrated  numerically  from  an  initially  subsonic 
condition  to  a  supersonic  one.  The  same  equations  and  notation,  except 
for  the  mass  flow,  given  in  Part  II  are  retained  here.  Thus,  the  quanti¬ 
ties  x,  y  and  a  of  Chapter  2.0  become  x  ,  distance  along  the  nozzle 
axis,  T  ,  temperature,  and  m  ,  mass  flow,  respectively.  Equation  (2-1) 
is  replaced  by 


^  P 
dx  "  Q  ' 

where 

Q  s  1  -  4  , 


and 


r  1  ^ 

Nj^  \  A  dx 
n, 

Ls  Pi  1 
i=l 


(5-1) 

(3-2) 


(3-3) 
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In  the  foregoing,  Mf  is  the  frozen  Mach  number, 

the  specific  heat  at  constant  pressure  of  species 


the  flow  speed, 
i 


n  the  mole  - 


ratio  of  species  i  ,  the  total  number  of  chemical  species,  and 


A  the  nozzle  area.  As  is  evident  from  definition  (3-2),  the  singularity 
occurs  at  Mp  =  1  .  At  this  point  P  must  also  be  zero  if  dT/dx  is  to 
have  a  finite  value.  Consequently,  the  singularity  occm’s  when  simulta¬ 
neously  Kf.  =  1  and  (f>  =  A"l(dA/dx)  .  The  quantity  (p  is  given  by 


N, 


I". 


1=1 


i=l 


j, 

^y.  ”i/ 

i=l  ' 


X 


Ll=l 


c  T  +  e  - 
p.  V.  D, 

■^1  1  i 


dn^ 

dx 


1  de 


i=l 


(3-4) 


where  R  is  the  universal  gas  constant,  the  vibrational  energy 

per  mole  of  species  1  ,  and  Sj.  the  energy  of  reaction  per  mole  of 

species  i  .  Equations  (3-1)  through  (3-4)  constitute  only  a  portion  of 
the  system  of  algebraic  and  ordinary  differential  equations  that  must  be 
solved  simultaneously  as  described  in  Part  II.  The  other  equations, 
however,  are  not  needed  here  since  only  equation  (5-1)  contains  a 
singularity . 


In  ordei’  to  clarify  the  subsequent  discussion,  the  nozzle  is  divided, 
as  shown  in  Figure  2,  into  five  regions.  Each  region  is  discussed  in 
turn  in  the  ensuing  sections.  Region  I,  discussed  in  Section  3.2,  con¬ 
sists  of  the  low-subsonic  portion  of  the  nozzle.  (Section  5.2  also  con¬ 
siders  how  an  initial  estimate  for  the  mass  flow  can  be  made.)  Region 
II  is  similar  to  region  I  except  that  XpQ  is  now  evaluated  at  each 
integration  step  in  order  to  locate  the  optimum  point.  Region  III  begins 
at  the  optimum  point  and  terminates  when  the  flow  becomes  supersonic. 

Since  this  region  contains  the  saddle  point,  an  inverted  form  of  equa¬ 
tion  (3-1)  is  used  here.  Region  IV  "patches"  the  altered  nozzle  area 
actually  used  in  the  forward  integration  In  region  III  to  the  given  nozzle 
area.  Region  V  Includes  the  remainder  of  the  nozzle.  The  final  section 
of  the  chapter.  Section  3.7,  discusses  how  the  mass  flow  is  modified  if 
a  more  precise  solution  is  desired. 

Three  input  parameters,  denoted  by  Vj^  (i  =  1,2,3)  ,  are  Introduced 
to  control  certain  arbitrary  features  of  the  optimum-point  method.  Thus 
X  =  Vp  initiates  region  II,  ^^  =  V2  initiates  region  IV,  and  Vg  con¬ 
trols  the  overall  accuracy.  The  quantity  V  defined  as  the  length  of 
region  IV  is  not  to  be  confused  with  the  input  parameters  Vp  . 
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3.P  REGION  I 

The  numerical  integration  of  the  noiiequilihriuin  nozzle-flow  eijua- 
tlons  starts  at  ,  which  for  simplicity  is  taken  as  zero.  At  this 
location,  the  flow  is  assumed  to  he  in  equlllhri'um  at  some  density 
and  temperature  T^  .  This  determines  the  initial  composition  n^^  ^ 
and  vlhratlonal  energy  e  .  As  in  Part  II,  the  nozzle  area  A’  is 

i>o 

considered  to  he  a  known  function  of  x  . 

The  mass  flow  m  and  stagnation  enthalpy  h^  ,  must  also  he  given. 
The  best  choice  for  m  ,  of  course,  would  he  the  critical  mass  flow  m„ 
Although  this  (ptantity  is  unknown,  it  satisfies  the  following  inetiual- 
Ities  (see  Bray  (1959)): 


m  <  m  < 
eq  cr  f 


where  nig^  is  the  equlllhrlum-flow  value  of  m  and  nif.  is  the  frozen- 
flow  value,  i.e.,  the  value  with  all  and  e^  held  constant  at 

their  initial  values.  Based  on  experience,  a  good  initial  guess  for  I, he 
mass  flow  is  0.95  nif  .  The  frozen  mass  flow  is  approximately  given  by 


1 


where  A^j^  is  the  nozzle  area  at  the  throat,  and  the  ratio  of  specific 
heats  is  given  by 


Equation  (3-5)  is  exact  when  the  initial  conditions  are  the  stagnation 

conditions.  Since  conditions  at  Xq=  0  correspond  to  a  non-zero  frozen 

Mach  number,  the  value  for  nij.  given  by  equation  (3-5)  is  somewhat  too 

small.  This  error  is  negligible,  however,  when  the  frozen  Mach  number 

at  X  =  0  is  small, 
o 
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Region  I  ends  when  x  =  V-j^  ,  where  the  Input  parameter  must  te 
chosen  sufficiently  large  that  the  first  minimum  of  i^pol  in  region  II 
corresponds  to  the  optltmitn  point,  fierierally  will  lie  at  a  nozzle 

location  slightly  upstream  of  the  throat. 

3.3  REGION  II 

After  each  integration  step,  XpQ  is  computed  and  stored  for  two 
steps.  Its  numerical  value  is  given  oy 


X 


PQ 


hjp+P'*'  Q  + 

?  i  +  ’  I  ! 
■|P-P  Q-QJ 


(3-7) 


where  h  is  the  integration  step  size  and  (  )  denotes  the  value  of 
(  )  at  x-h  .  Equation  (3-7)  is  thus  a  finite -difference  form  of 
equation  (P-5)*  A  minimum  of  l^pql  occurs  when  both 


I^pqI  s  |Xp|^|  , 


(3-8) 


are  satisfied,  where  XpJ  is  the  value  of  XpQ  at  x-  2h  .  Conditions 
(3-8)  are  checked  after  each  integration  step  in  region  II.  When  satis¬ 
fied,  the  value  of  x  corresponding  to  is  the  optimum  point  . 

At  this  time,  the  actual  integration  is  at  Xqp^+  h  .  Since  region  III 
starts  at  Xgp^.  ,  it  is  advisable  that  the  quantities 

T,  dT/'dx  ,  ti^(i=l,...,Nj^)  ,  e^  (i=.1 ,  ,N^)  , 

be  si,ored  for  one  step  in  region  II. 

3.)|  REGION  III 

Region  III  begins  where  x  =  and  continues  until  the  frozen 

Mach  number  is  greater  than  unity.  The  temi)crature  Tjyj  in  this  region 
is  given  by  (see  equation  (2-6)) 


(3-9) 


As  pointed  out  in  Chapter  2.0,  it  is  desirable  to  retain  equation  (3-1) 
in  region  III.  This  is  accomplished  by  inverting  this  equation  to  obtain 


dA 


III 


dx 


c  n. 
Pi  1 


(3-10) 
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r< 

ii 


1 ", 
:  •* 

ii 


In  Lhis  eijuation  Ajj-p  is  an  unknown  noazlc  area.  In  region  III  (and 
later  in  region  IV)  the  given  nozzle  area  is  denoted  by  Ag  to  distin¬ 
guish  it  from  the  nozzle  area  actually  used  in  the  forward  integration. 
These  latter  nozzle  areas  are  denoted  by  Ajjj  and  Ajy  .  In  the  above 
form,  equation  (5-10)  contains  no  singularity. 

Equations  (3-9)  aad  (3~10)  in  conjunction  with  the  other  nozzle-flow 
equations  yield  an  exact  numerical  solution  providing  A™  is  replaced 
by  Ajjj  and  Ajy  in  regions  III  and  IV.  In  other  words ^  the  mass  flow 
m  being  used  becomes  the  critical  mass  flow  for  a  nozzle  whose  area 
distribution  is  the  given  (xie  in  regions  I  and  II  and  is  specified  by 
Ajii  and  Ajy  in  regions  III  and  IV.  A  convenient  measure  of  the  accu¬ 
racy  of  the  optimum-point  method  that  is  useful  for  mass -flow  modifica¬ 
tion  discussed  in  Section  3.7,  is  thus  given  by  the  percent  error 

A  -  A__^  „ 

E  s  ^  ^  X  10^  ,  (3-11) 


Region  III  Is  terminated  when  ^  ^  where  the  input  parameter 

V^  is  slightly  greater  than  unity  (supersonic  flow).  This  nozzle  loca¬ 
tion  is  referred  to  as  the  transfer  point  . 

3.5  REGION  IV 

At  the  transfer  point,  the  area  Ajjj  and  its  derivative  dAjjj/dx 
do  not  equal  Ag  and  dA  /dx  .  Consequently,  region  IV  constitutes  a 
short  nozzle  se^ient 

X.  S  X  S  X,  +  V  (3-12) 

tr  tr  ' 

that  is  used  to  "patch"  Aqjq  to  the  given  area  Ag  at  x^j,+  V  . 

How  the  positive  quantity  V  is  found  is  described  near  the  end  of  this 
section.  Thus,  Ajy(x)  is  used  to  join  Ajjj  gradually  to  A„  .  The 
nonequilibrium  equations  of  regions  I  and  II  are  used  here.  This  section 
therefore  describes  a  procedure  for  calculating  Ajy  . 

Considerably  more  effort  was  expended  in  this  phase  than  in  the  rest 
of  this  work  combined.  A  satisfactory  but  sophisticated  procedure  for 
finding  Apy  was  arrived  at  only  after  an  exhaustive  search.  Any  simpler 
procedure  may  work  for  specific  stagnation  conditions  and  given  nozzle 
geometry  but  will  prove  inadequate  for  others.  Thus,  a  wide  variety  of 
stagnation  conditions  and  given  nozzle  geometries  were  used  to  check  the 
various  procedures  tried.  Only  the  procedure  given  here  proved  to  be  sat¬ 
isfactory  in  ail  cases.  For  these  reasons,  the  nature  of  the  difficulty 
will  first  be  described  in  detail. 
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At  the  transfer  point,  the  qiuantity  Ajjj(dAjj-|-/dx)  -<p  in  equation 
(3-3)  is  close  to  zero  since  it  is  zero  just  upstream  of  this  point. 
Inasmuch  as  Ajy(dAjy/dx)  and  -  <j>  are  nearly  equal  at  the  start  of 
region  IV,  a  small  error  in  the  value  of  A£^(dAjy/dx)  will  result  in 
a  large  error  in  dl'/dx  .  Therefore,  unless  Ajy  matches  Ajjj  at 
in  a  smooth  fashion,  the  temperature  gradient  in  region  IV  will 
differ  drastically  from  that  in  region  III,  i.e.,  from  (dT/dx)Qp^.  . 


To  avoid  this  discontinuity  the  three  quantities 


1  ^  ^/l 

A  dx  ^  dx\A  dx  ..’  •’ 


(3-13) 


are  required  to  be  continuous  at  x-i^p  .  These  conditions  are  generally 
sufficient  to  assure  a  smooth  variation  of  T  and  dT/dx  at  x^^  . 

The  area  function  A“^(dA/dx)  is  the  most  important  of  the  fore¬ 
going  quantities  since  it  appears  directly  in  equation  (3-1).  Its 
variation  with  x  in  the  vicinity  of  the  singularity  is  sketched  in 
Figure  3.  In  this  figure  the  solid  curve  represents  the  given  area 
function  Ag  (dAg/dx)  ,  while  the  dashed  curve  labeled  III  represents 
Ajji(dAiii/ax)  .  The  dashed  curve  corresponds  to  a  mass  flow  somewhat 
greater  than  the  critical  one  since  dAj-jj/dx  is  greater  than  dAg/dx  . 
If  the  mass  flow  is  less  than  the  critical  one  the  dashed  curve  would  be 
below  the  solid  one.  The  other  curves  in  Figure  3  are  described  later. 


The  continuity  requirements  (3-13)  do  not  always  completely  remove 
the  discontinuity  in  dT/dx  at  X|.j.  .  When  Runge-Kutta  integration  is 
used  (see  Part  III  for  further  details),  an  unstable  feedback  mechanism 
may  occur  during  the  first  integration  step  in  region  IV.  This  mecha¬ 
nism  operates  as  follows:  As  already  noted,  a  small  error  in  the  value 
of  Ajy(dAj-y/dx)  first  results  in  a  larger  error  in  dT/dx  .  Subse¬ 
quently,  the  value  of  the  temperature  will  be  in  error  and  thereby  cause 
errors  in  the  values  of  dn^/dx  and  de„  /dx  .  The  value  of  0  ,  given 

i 

by  equation  (3-4),  therefore  will  also  be  somewhat  in  error.  If  the 
errors  In  the  values  of  Ajy(dAjy/dx)  and  0  are  additive,  as  is  occa¬ 
sionally  the  case,  then  the  resulting  value  for  dT/dx  will  now  be 
appreciably  in  error.  This  leads  to  further  errors  in  the  values  of 


dn^/dx 


etc . 


The  above  mechanism  operates  only  during  the  first  Runge-Kutta 
integration  step  in  region  IV,  no  matter  how  small  the  size  of  the  inte¬ 
gration  step.  Basically  these  difficulties  are  due  to  the  nozzle  area 
not  being  analytic  at  x^p  .  Matching  third  or  higher  derivatives  of 
the  area  at  x^^  woi,ild  partially  alleviate  the  problem.  This  is  im¬ 
practical,  however,  as  will  become  apparent  later.  The  discontinuity 
in  dT/'dx  ,  which  is  minimized  when  the  procedure  described  later  is 
used,  decays  rapidly  after  the  initial  integration  step  in  region  IV. 
Thus,  by  the  time  three  integration  steps  are  completed  the  temperature 
gradient  has  returned  to  a  value  consistent  with  (dT/dx)Qp^  . 
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The  feedback  difficulty  just  described  is  not  alleviated  by  increas¬ 
ing  V?  .  Increasing  Vg  merely  shifts  the  transfer  point  down¬ 

stream  and  allows  the  error  E  to  increase  further. 


It  is  instructive  to  describe  what  happens  when  a  polynomial  is 


used  to  represent  A- 


IV 


,  since  this  formulation  is  quite  simple.  For 
lOlynomial  in  x  could  be  used  to  match  A 
at  and  the  same  quantities  at  x^j,~+V 


III 


Al^ ( dA-ry/dx )  is  sketched  in  Figure  3  as 


example^  a  fifth-degree  n 
dAiii/dx  and  d^Ajjj/dx^ 

The  resulting  variation  of  ^  ll-A.  y  X  O  IVt  U  1  It;  VX  Xll  IX^ 

the  dotted  curve  labeled  "polynomial  area  variation."  This  curve  con¬ 
stitutes  a  poor  approximation  to  the  given  one  (solid  curve)  and  will 
result  in  values  for  dT/dx  considerably  different  from  (dT/dx)Qp^ 
throughout  most  of  region  IV. 


As  mentioned  earlier,  when  the  mass  flow  is  greater  than  the  crit¬ 
ical  one,  Ajjj  and  dA^^^/dx  at  x^^,  will  generally  be  greater  than 
their  corresponding  given  values.  Thus,  if  Ajy  is  to  match  A_  at 
xtr  tV  ,  then  the  gradient  dAjy/dx  must  be  less  than  dAg/dx  curing 
part  of  region  IV.  This  is  the  reason  why  any  area  function 
Ajy(dAjy/dx)  must  cross  Ag^(dA  /dx)  (see  Figures  3  and  5).  Similar 
reasoning  also  leads  to  the  conclusion  that  the  two  curves  cross  when 
the  mass  flow  is  less  than  the  critical  one. 


We  shall  now  describe  the  method  of  formulation  and  the  procedure 
finally  arrived  at  to  determine  Ajy  .  The  conditions  Imposed  are  as 
follows: 

(a)  At  X  =  x^j,  ,  we  require  that 


\v  “  ■^III 


(3-14) 


dA, 


IV 


dx 


1 

A 

III 


dA 


III 


dx 


(3-15) 


dx  yA^y  dx 


_d_  1 

dx  \  Ajjj 


(3-16) 


(b)  At  x  =  x^p  +  V  ,  we  require  that 


"iv  -  \  ' 


(3-17) 


A^,,  dx  A  dx 
IV  g 


(3-18) 
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,  /  ,  dA^A  j  /  1  dA  \ 

_1 _ IV 1  =  i_i  1 _ s\ 

^\^IV  J  ^  I  ‘ 


(3-19) 


(c)  In  region  IVj  we  also  reiiulre  that  Aji(dAjY/dx)  be  a  reason¬ 
able  approximation  to  A"^(dA  /dx)  .  This  condition  will  be 
formulated  more  precisely  later. 


In  addition  to  x  ,  the  variable  x  is  also  used,  where 


X 


=  y(x-Xtr)  -1  . 


(3-20) 


Equation  (3-20)  linearly  transforms  the  interval  (x^j,,x^^+V)  into 
(-1,1)  .  Since  A£y(dA^y/dx)  ,  rather  than  A^^  ,  Is  of  primary  Impor¬ 
tance,  the  following  formulation  is  used: 


dA 


IV 


dx 


dA 


(3-21) 


where  and  Ag  are  considered  to  be  functions  of  x  and  the 

functlon'^^\^  Is  still  to  be  found.  To  determine  Ajy  ,  equation  (3-21) 
Is  integrated  from  x^^.  to  x  ,  thereby  resulting  in 


^iv(^) 


^  A^(X)  , 

g'  tr' 


(3-22) 


where 


^^(x)  s  IJ  K^(x)di( 


(3-23) 


Condition  (3-lh)  is  then  satisfied  by  replacing  Ajy(x^y)  in  equation 
C3-22)  by  the  known  quantity  Ajjj(x^j.)  .  Condition  (3-15)  becomes 
\^(-l)  =  ,  where  the  known  constant  is  given  by  (see  Figure  3) 


x=x 


tr 


(3-24) 
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15  1?  „  . 

■i3= 


x=x 


^  ’ll  ^  3  ’Ig 


tr 


(3-31 ) 


The  foregoing  results  satisfy  conditions  (a)  and  (b)  given  by  equa¬ 
tions  (3-ll|)  through  (3-19).  Condition  (c)  is  approximately  satisfied 
by  choosing  fg  appropriately.  For  simplicity,  it  is  taken  as 


f3(x)  =  TIjj^X  , 


(3-32) 


where  the  constant  is  still  to  be  determined. 

Figure  3  shows  a  sketch  designated  by  (iv)  of  Ajy(dA]-y/dx)  for 
a  positive  value  of  ri^  ,  while  Figure  ^l-  shows  the  corresponding 
ki/hi  curve.  Figure  5  is  similar  to  Figure  3  except  that  tj-j^  is 
negative.  In  either  case,  Xp  has  an  extremum  in  the  vicinity  of  the 
origin.  This  extremum  and  its  location  are  designated  by  Xp  and  x^  . 
As  is  evident  from  Figures  3  and  4,  condition  (c)  is  approximately  ful- 
I’llled  if  |Xp|  is  minimized. 

For  fixed  values  of  hp,  tl2  >  and  qg  ,  x*  is  determined  by  the 
condition 

SX-, 

=  0  ,  (3-33a) 

die 


or 


■11^  =  f  -ip  --  i 

To  minimize  |xT|  ,  we  first  note  that  Xp  =  Xp(x*(iij^)  jpjj)  .  Hence, 
tlie  desired  minimum  occurs  when 


'>'11,  ’Si/. 

x=x* 


=  0  . 


(3-34a) 


The  second  term  on  the  right-hand  side,  however,  is  zero  by  virtue  of 
equation  (3-33a).  Thus,  equation  (3-34a)  can  he  shown  to  reduce  to 

(x*+l)^(x*  -l)^x*  =  0  .  (3-34b) 


The  solution  of  equations  (3-33b)  and  (3-34b)  readily  results  in 
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X*  =  0  , 


% 


3 


(3-35) 


When  equations  (3-28)  through  (3-32) 
hined,  the  following  Is  obtained; 


and  equations  (3-35)  are  com- 


(3-36) 


The  parametei'  V  is  now  determined  by  setting  -(\^/riq)  =  (3/^)  in 

equation  (3-36)  where,  for  convenience,  the  term  1^2/(161^2^)  is  neglected. 
When  examined  a  posteriori,  the  value  of  this  term  is  indeed  found  to  be 
small  compared  with  unity.  Thus,  V  is  given  by 


V 


6_ 

’ll 


&n 


tr 


(3-37) 


A  smaller  value  than  3/4  for  ■(4*/'')l)  would  result  in  a  larger  value 
for  V  ,  while  a  larger  value  than  3/4  would  result  in  a  smaller  V  . 
The  choice  made  here,  however,  proved  to  be  an  acceptable  one.  As  a  con¬ 
sequence  of  equation  (3-37),  the  constant  tig  is  given  by 


i  \  • 


(3-38) 


The  pertinent  equations  of  this  section  are  summarized  in  a  form 
suitable  for  numerical  work  in  Appendix  I. 

3.6  REGION  V 

This  region  includes  the  remainder  of  the  nozzle.  The  given  nozzle 
area  is  used,  and  the  fojrward  integration  proceeds  as  described  in 
Part  II. 

3.7  MASS-FLOW  MODIFICATION 

The  initial  value  for  the  mass  flow  when  chosen  according  to  the 
method  in  Section  3.2  may  differ  by  as  much  as  5^^  from  the  critical  mass 
flow.  To  obtain  a  more  precise  calculation,  a  method  is  now  described 
for  modifying  the  mass  flow  and  the  stagnation  enthalpy. 

After  each  integration  step  in  region  III  the  following  check  Is 
performed: 


|E|  ^  Vg  , 


(3-39) 
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where  E  is  defined  by  equation  (3-11)  and  V3  is  a  positive  input  param¬ 
eter  that  controls  the  accuracy  of  the  method.  As  long  as  condition  (3-39) 
is  satisfied,  the  integration  proceeds  as  described  in  Section  3.1)-.  If 
condition  (3-39)  is  not  satisfied,  then  the  integration  starts  over  again 
at  0  with  new  values  for  the  mass  flow  and  stagnation  enthalpy. 

To  determine  a  new  value  for  the  mass  flow  that  is  closer  than  the 
preceding  one  to  the  critical  mass  flow,  two  different  fictitious  frozen 
flows  are  compared.  The  first  flow  uses  the  current  mass  flow  1%  ,  the 
temperature  and  nozzle  area  at  0  ,  and  the  temperature  and  nozzle 

area  at  the  point  in  region  III  where  (3-39)  is  not  satisfied.  Accord¬ 
ing  to  the  usual  frozen-flow  equations,  is  given  by 


1, 

2 


III 

T 


m  =  c  • 
a  o . 


(3-40 ) 


:  '  A  T 
i  O  ,  I  o  ' 

'7iiV 


^0-1 


-  1 


where  is  defined  by  equation  (3-6)  and  the  constant  Cp  is 


1 
.  2 


wv 


= 


y  \  v'l  I 

-7;  RT  )  n  ;  p  A  . 

-J-/  o  /Li  1^0  '^00 

i=l  .! 


(3-41) 


The  second  flow  is  similar  to  the  first  except  that  the  improved  muss 
flow  replaces  m^  and  the  value  of  Ag  at  the  point  where  (3-39) 

is  not  satisfied  replaces  Ajjj  .  Thus  is  given  by 


%  = 


J — 

iLm  -| 

0 

2 

k  \ 

0 

1 

0 

1  °  1  1 

iWii/ 

(3-42) 


Taking  the  ratio  of  equations  (3-4o)  and  (3-42)  then  results  in  the 
desired  mass -flow  correction  formula 
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1. 

2 

I 


■ 

■,  ""o  I 


'  T 


Lv^’iii 


A 's'' 

_£1 

V 


A  \ 

'^iii; 


(3-43) 


The  stagnation  enthalpy  h.(-  is  also  modified  in  order  that  the 
new  nozzle -flow  solution  can  start  at  the  same  equilibrium  state 
previously  used.  Thus,  p^,  T^,  n^  ^  ,  and  e.^  remain  unchanged. 

^  i  ,  o 

This  is  readily  accomplished  providing  is  unchanged,  since  the 

other  quantities  are  initial  conditions.  The  requirement 


*^o  A  (u  ) 
o'  o' c 


(U  ) 

o'. 


(3-44) 


is  therefore  imposed,  where 


(u_ )  is  the  flow  speed  at 
'  °  a 


x^=  0  when 


the  maos  flow  is  m  and  (u  )  is  the  corresponding  flow  speed  when 

b 

the  mass  flow  is  •  The  energy  equation  is  now  evaluated  at  0 

for  flows  (a)  and  (b)  as  follows: 


h 


o 


(31: 5a) 


vh),  -o  *  I  (v^  ' 

where  h^  is  the  enthalpy  at  Xq=  0  .  By  combination  of  equations 
(3-1(4)  and  (3-45),  the  desired  stagnation-enthalpy  correction  formula 
Is  obtained  as 
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4.0  RESULTS  AND  DISCUSSION 
4,1  GIVEN  CONDITIONS 

This  chapter  describes  the  results  of  one  noneijuilibrium  nozzle- 
flow  calculation  that  utilized  the  foregoing  optimum-point  method. 

This  calculation  and  all  others  that  were  performed  are  based  on  the 
same  gas  model  used  in  Parts  I  and  II.  This  is  a  model  of  air  consist¬ 
ing  of  the  five  species  0,  N,  NO,  N2  ,  and  Og  (identified  by  sub¬ 
scripts  i  =  1,...,5  ,  respectively).  The  same  eight  reactions  and  the 
same  rate  constants  are  also  retained. 

The  variation  of  nozzle  area  with  distance  for  all  five  regions  is 
given  by 

A  =  10  -  6x  +  .  (4-1) 

g 

This  equation  represents  a  hyperbolic  nozzle  with  a  throat  area  of 
1  cm^  located  at  x  =  3.0  cm  .  All  numerical  integrations  started  at 
X  =  0  cm  where  the  area  is  10  cm^  and  concluded  at  x  =  4  cm  where 
the  supersonic  flow  is  well  into  region  V. 

At  the  Initial  point  x  =  0  cm  the  flow  is  assumed  to  be  in  equl 
librium  at  a  temperature  of  8000°K  and  a  density  of  2.8o  x  lO"^  gm/cm 
These  conditions  were  chosen  so  that  a  significant  portion  of  the  dia¬ 
tomic  species  would  be  dissociated  at  x  =  0  .  In  addition,  to  insure 
that  the  flow  would  depart  from  equilibrium  upstream  of  the  throat,  a 
small  initial  density  as  compared  with  that  used  in  Parts  I  and  II  was 
necessary.  For  all  practical  purposes,  the  foregoing  are  also  the 
stagnation  conditions  since  the  initial  frozen  Mach  number  Is  0.0^48. 
The  equilibrium  composition  and  vibrational  energies  corresponding  to 
this  density  and  temperature  are  given  in  the  following  table: 


EQUILIBRIUM  COMPOSITION  AND  VIBRATIONAL  ENERGIES 
FOR  AIR  AT  8OOOOK  AND  2.8o  X  10 "5  GM/CM^ 


1 

SPECIES 

n^(M0LES/GM) 

ey^(DYNB  CM/MOLE) 

1 

0 

.14598473x10'^ 

2 

N 

. 45508697 X 10'^ 

- 

3 

NO 

.60071 580 X  lO'^ 

. 55920669  X 10^^ 

4 

No 

.45156159x10'^ 

.53593265x10^^ 

5 

.72687672  X 10'^ 

.5"'674503x  10^^ 
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initial  values  for  the  mass  flow  and  sLagnation  enthalpy  were 
chosen  in  accord  with  the  discussion  in  Chapter  3,  The  final  values  for 
these  quantities  are  given  in  the  next  section.  The  following  values 
were  used  for  the  input  parameters  : 

=  2.8  , 

Vg  =  1.01  , 

^3=  • 


The  optimum- point  method  was  also  used  in  conjunction  with  the 
approximate  equilibrium  formulation  given  in  Part  II  and  modified  in 
Appendix  II  of  this  report.  The  same  nozzle  geometry  and  initial  con¬ 
ditions  as  described  above  were  used  (except  for  mass  flow  and  stagna¬ 
tion  enthalpy).  The  only  changes  were  that  the  nj^  were  sequenced  as 
follows;  O2,  N,  NO,  Ng,  0  .  The  reason  for  this  change  is  discussed 
In  Appendix  II.  In  addition,  for  reasons  given  in  the  next  sectirn,  a 
value  of  1.0  was  used  for  . 

The  calculations  presented  here  were  chosen  because  they  proved  to 
be  the  most  troublesome  in  terms  of  the  instability  described  in  Section 
3.5.  In  fact,  of  the  many  procedures  for  determining  Ajy  that  were 
tried  only  the  one  given  here  dealt  successfully  with  this  case.  In 
terms  of  instability,  the  approximate  equillhrlura  caDculation  was  more 
troublesome  than  the  nonequilibrium  one.  This  point  is  discussed  further 
in  Section  k.y. 

4.2  RESULTS  01'  THE  OPTIMUM- POINT  METHOD 

On  the  basis  of  equation  (3-5),  the  foregoing  conditions  result, 
in  a  frozen  mass  flow  of  4.1925  gm/sec.  Although  equation  (3.5)  is 
approximate,  this  value  in  very  close  to  the  true  value  since  the  initial 
Mach  number  is  nearly  zero.  The  Initial  value  for  the  nonequilibrium 
mans  flow  was  Uiex-efore  chosen  to  be  0.95x  4.1925  =  3.98  gm/sec.  The 
final  values  for  the  mass  flow  and  stagnation  enthalpy  came  out  to  be 
4.0540000  gm/sec  and  7-00648?? x 10^®  cm^/sec^  ,  respectively.  Figures 
5  thi’ough  11  are  based  on  these  nonequilibrium  values.  This  value  for 
the  mass  flow  Is  slightly  helow  the  critical  one,  as  is  evident  from 
Figure  5-  A  subsequent  calculation  with  a  mass  flow  of  4.0830000  gm/sec 
showed  that  this  larger  value  was,  in  comparison,  well  above  the 
critical  one.  Thu.s,  the  mass-flow  value  of  4,054  gm/sec  is  accurate  at 
least  to  three  significant  figures. 

With  regard  to  accuracy,  let  us  first  compare  the  location  and 
area  of  the  throat  for  the  two  nozzle  areas,  Ag  and  .  This  com¬ 

parison  is  presented  in  the  following  table: 
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LOCATION  AND  AREA  OF  THE  NOZZLE  THROAT 


NOZZLE  THROAT  THROAT 

AREA  LOCATION  AREA 


A  3.0000  cm  1.00000  cm^ 

g 

Ajjj  3.0083  cm  0.9993?  cm^ 


Another  pertinent  ijuantltjr  is  the  maximum  value  of  the  percent  error  E  . 
This  occurs  at  and  is  0.21?^,  well  below  its  upper  limit  of 

0.5^.  Finally,  we  note  that  the  calculated  area  function  A“^(dA/dx)  , 
shown  in  Figure  5,  differs  but  slightly  from  the  corresponding  given 
area  function  in  both  regions  III  and  IV.  (See  Figure  3  for  a  more 
detailed  description  of  the  curves  in  Figure  5- )  We  therefore  conclude 
that  Ajjj  and  A^^  do  not  differ  significantly  from  Ag  for  reason¬ 
ably  accurate  values  of  the  mass  flow,  i.e.,  three  significant  figures. 
Calculations  based  on  a  more  accurate  estimate  for  the  mass  flow  than 
4.05^  gm/sec  could  be  readily  obtained.  Such  a  calculation,  hovjever,  was 
not  deemed  necessary  since  the  results  for  4.09^  gm/sec  and  4.o83  gm/sec 
differed  only  negligibly.  Thus,  when  the  optimum-point  method  is  used 
the  resulting  solution  differs  negligibly  from  the  critical  one  for 
reasonably  acciirate  values  of  the  mass  flow. 

For  the  approximate  equilibrium  calculation,  the  values  for  mass 
flow  and  stagnation  enthalpy  are  3.7268000  gm/sec  and  7.0048583x  10^*^ 
cm^/sec^,  respectively.  This  calculation,  which  is  also  shown  in 
Figures  6  through  11,  differed  from  the  nonequilibrium  one  in  two  respects. 
First,  greater  accuracy  was  required  for  the  mass-flow  value.  This  is 
basically  a  consequence  of  the  curves  for  P  =  0  and  Q  =  0  being  close 
together  upstream  of  the  saddle  point.  As  a  result,  for  any  value  of  the 
mass  flow  different  from  the  critical  one,  the  calculation  behaves  as  if 
the  singularity  were  located  at  the  geometric  throat,  where  the  equilib¬ 
rium  Mach  number  is  unity,  rather  than  at  =  1  ,  The  extent  of  region 

III  is  therefore  somewhat  longer  and  the  percent  error  at  the  transfer 

point  is  larger  as  compared  with  the  nonequilibrium  calculation.  Hence, 

a  larger  value  was  used  for  Vj  in  the  equilibrium  calculation.  As  in 

the  nonequilibrium  case,  however,  the  solution  for  slightly  different 
values  of  the  mass  flow  (i.e.,  m  =  3.7271  gm/sec  )  differed  negligibly 

from  the  one  given  here.  Thus,  the  mass-flow  value  of  5.7268  gm/sec  is 

accurate  to  at  least  four  significant  figures. 

The  second  difference  is  that  the  instability  at  the  start  of  region 

IV  is  more  pronounced  for  the  equilibrium  calculation.  This  resulted 
in  a  reduction  in  the  integration  step  siiie  at  the  start  of  region  IV. 
Since  the  Instability  decays  rapidly,  however,  the  step  size  then  proceeds 
to  increase  to  its  former  level. 
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h . 3  GASDYWAMIC  KESULTS 

Figures  6  through  11  show  the  results  of  the  nonequilibrium  and 
equilibrium  calculations.  In  all  the  figures  the  abscissa  is  the  dis¬ 
tance  along  the  nozzle  axis. 

'Figure  6  shows  the  temperature  variation.  For  Interest  the  temper¬ 
ature  variation  is  also  shown  for  a  flow  frozen  at  the  same  initial 
conditions.  From  this  figure  it  is  evident  that  the  nonequilibrium 
solution  starts  to  diverge  from  the  equilibrium  one  at  about  x  =  1.8  cm 
where  the  frozen  Mach  number  is  0.23h.  Along  the  entire  nozzle ^  the 
nonequilibrium  and  frozen-flow  temperatures  are  quite  close,  indicating 
that  once  the  flow  diverges  from  the  equilibrium  one  it  freezes  rapidly. 
This  result  is  not  a  general  one  but  rather  is  a  consequence  of  the 
specific  stagnation  conditions  and  nozzle  geometry  used  for  these 
calculations.  In  particular,  a  rapidly  varying  nozzle  geometry,  which 
is  the  case  here,  tends  to  favor  rapid  freezing.  The  large  difference 
between  the  equilibrium  and  nonequilibrium  temperature  after  x  =  cm 
is  due  to  nitrogen  recombination  in  the  equilibrium  flow.  This  recom¬ 
bination  is  briefly  discussed  later. 

Figure  7  shows  the  density  and  pressure  variations.  In  comparison 
with  Figure  6,  it  is  seen  that  density  and  pressure  are  not  as  sensitive 
as  the  temperature  to  nonequilibrium  effects.  This  result  is  in  accord 
with  other  nozzle-flow  calculations,  such  as  those  of  Bray  (l959). 

Figure  8  shows  the  variation  of  the  flow-speed.  The  equilibrium 
and  nonequilibrium  speeds  differ  slightly  at  x  =  0  cm  because  the 
respective  mass  flows  are  different.  Figures  9  tail's  10  show  the  varia¬ 
tion  of  chemical  composition  and  vibrational  energy,  respectively. 

All  quantities  shown  in  Figures  6  through  iO,  with  the  exci-.'pt.lon 
of  flow-speed,  are  nearly  constant  in  the  I'egion  from  x  ^  0  cm  to 
X  =  1.8  cm  .  In  this  low-subsonic  region,  the  nonequilibrium  (or 
equilibrium)  solution  is  given  approximately  by 

uA  =  ™  =  constant  ,  (4-2) 


where  all  other  variables  ai'c  constant.  Thus,  the  nonequilibrium  and 
equilibrium  solutions  will  be  close  here  regardless  of  the  values  of 
the  chemical  and  vibrational  relaxation  lengths.  In  this  region,  all. 
relaxation  lengths  are  directly  proportional  to  velocity  and  therefore 
increase  steadily.  Should  any  of  the  charai^teristic  lengths  become 
sufficiently  large  here,  the  flow  will  Lhen  start  to  diverge  from  equi¬ 
librium  as  soon  as  compressibility  becomes  Important.  This  occurred 
in  the  calculation  shown  here. 
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Approximate  values  for  the  characteristic  lengths  at  x  =  0  cm  and 
at  X  =  1.8  cm  were  calculated  according  to  the  theory  In  Part  II  and 
In  Appendix  II.  They  are  presented  In  the  following  table: 


CHEMICAL  AND  VIBRATIONAL  CHARACTERISTIC  LENGTHS 


RATE  EQUATION 

FOR 

CHARACTERISTIC 
LENGTH,  CM, 

AT  X  =  0  cm 

CHARACTERISTIC 
LENGTH,  CM, 
AT  X  =  1.8  cm 

0  , 

CHEMICAL 

1.32  XlO"'^ 

5.4o  xlO"^ 

N  , 

tt 

6.71x10'^ 

2.1k  xlO“^ 

NO  , 

It 

-h 

5.07x10 

-3 

2.07X  lO"^ 

-3 

NO  , 

VIBRATIONAL 

1.30X10 

5-32  xio 

No  , 

tt 

1.59X10"^ 

6.50x10"^ 

4. 

tt 

1.11  xio'^ 

4.55x10"^ 

At  X  =  0  cm  all  of  the  characteristic  lengths  are  small.  At 
X  =  1.8  cm  ,  however,  the  mole-mass  ratio  for  atomic  nitrogen  should 
start  to  diverge  from  its  esiuilibrium  value,  whereas  the  mole-mass 
ratios  for  atomic  oxygen  and  nitric  oxide  are  still  close  to  their 
equilibrium  values.  Since  the  mole-mass  ratios  for  diatomic  nitrogen 
and  diatomic  oxygen  are  given  by  the  equations  for  conservation  of 
components,  the  former  will  also  start  to  diverge  from  equilibrium 
while  the  latter  will  not. 

The  preceding  conclusions  are  not  discernible  from  Figure  9- 
Indeed,  this  figure  indicates  that  nitric  oxide  is  the  first  species  to 
diverge  from  equilibrium  and  that  atomic  nitrogen  remains  close  to  its 
equilibrium  value  until  x  =  S.O  cm  .  These  fallacious  conclusions  are 
a  direct  I’esult  of  the  log  scale  used  for  the  chemical  composition,  a 
convi’nlout  and  common  practice.  To  clarify  the  situation,  the  variation 
of  1  n gjjl  is  shown  in  Figure  11,  where  is  tho  mole-mass 

ratio  of  ’species  1  for  the  equilibrium-flow  solution.  This  figure 
clearly  shows  that  the  mole-mass  ratios  for  atomic  and  diatomic  nitrogen 
diverge  from  their  equilibrium-flow  value  far  more  appreciably  than  do 
those  for  the  other  species.  Figure  11  is  thus  in  accord  with  the 
findings  based  on  the  characteristic  lengths.  The  valley  in  the  curve 
for  diatomic  oxygen  in  Figure  11  is  due  to  the  equilibrium  and  nonequi¬ 
librium  mole-mass  ratios  crossing,  as  shown  in  Figure  9* 

According  to  the  foregoing  table  of  characteristic  lengths,  the 
vibrational  energy  of  nitric  oxide  will  diverge  from  its  local  equilib¬ 
rium  value  upstream  of  the  other  vibrationai  energies.  Thus  Figure  10, 
while  correct,  cannot  be  used  to  determine  now  close  any  vibrational 
energy  is  to  its  local-equilibrium  value.  This  situation  has  been 
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\i 


discussed  In  detail  in  Part  II. 

Figure  9  shows  that  a  significant  amount  of  nitrogen  recombination 
occurs  in  the  equilibrium-flow  solution.  This  recombination  accounts 
for  the  large  difference  between  the  equilibrium  and  nonequilibrium 
temperatures.  In  the  equilibrium  calculation,  the  mole-mass  ratios  for 
the  other  species,  i.e.,  0,  Og  ,  and  NO  ,  remain  fairly  constant. 

This  is  due  to  the  small  decrease  in  the  equilibrium  temperature,  which 
favors  recombination,  being  counterbalanced  by  the  large  decrease  in 
the  equilibrium  density,  which  favors  dissociation. 
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APPKWDTX  I 

ANALYTICAL  EXPKESSIONS  FOR  Ajy(x) 


The  variation  of  nozzle  area  with  dictancc  in  region  IV  is  given 
helow.  At  the  start  of  this  region  the  parameters  tjj,  V  ,  and  ijg 
are  calculated.  Equations  (l-5)  through  (l-9)  then  determine  A_„ 
and  dAjy/dx  . 


>^1  - 


:  _L  ^ 

•  A  .J.J  dx  '  A  dx 

^  -  x=x,  -h 

tr 


(i-1) 


'  1  ^^III  ^  ^ 

A^II  dx  -  Ag  dx 
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tr 
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1  -  =^=^tr 
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(i-M 
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\^(x)  =  4(3x^i-  x'^-  6x  -3)t)^  +  (x  +l)(4x®-  X  -l)Tipl  (l-6) 
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APPENDIX  II 

MODIFICATION  OF  THE  APPROXIMATE  EQUILIBRIUM  METHOD  OF  PART  II 


Two  modifications  of  the  approximate  eiiuilibrium  method  of  Part  II 
are  given  here.  The  need  for  these  changes  became  apparent  when  equi- 
librliun  nozzle-flow  solutions  were  sought  that  required  the  optimum- 
point  method.  Once  these  changes  were  Incorporated,  along  with  the  con¬ 
dition  that  Pii-%j[  /  0  (l  =  1,...,N^)  as  explained  in  Part  II,  no 
further  difficulties  were  encountered.  (See  Part  II  for  definition  of 
symbols. ) 

The  first  change  requires  that  y  be  replaced  by  -  7^  in  equation 
( SLc )  of  Part  II  as  follows : 


dx 

where  the  |7^ |  are  large  constants.  To  understand  the  reason  for  this 
change,  we  first  derive  the  characteristic  lengths  for  equations 

(ll-l).  The  derivation  Is  explained  In  detail  in  Appendix  III  of  Part  II. 
Thus,  the  are  given  by 


r  i 


i  — 


(ii-i) 


-1 


i  =  1,...,N, 


(11-2) 


with  density,  temperature,  and  all  nj  except  n^  held  fixed.  The 
quantity  F^  is  defined  by 


F 


i 


'.n. 

1 


i=l,...,Ng.  (II-3) 


When  the  characteristic  lengths  are  defined  as  above  (i.e.,  without 

absolute  value  signs  and  with  density  held  fixed),  they  are  equivalent 
to  given  in  Part  III.  After  the  partial  differentiation  is  per¬ 

formed,  and  conservation  of  components  as  well  as  Lj^  ~  1  (.law  of  mass 
action)  are  accounted  for,  we  obtain 


e 


ci 


1 


i=l,...,N3.  (II-4) 


29 


AREA 


Fig.  3  Sketch  of  the  Area  Function  A  ^  (dA/dx) 
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Fig.  6  Temperature 
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Fig.  7  Pressure  and  Density 
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Fig.  8  Speed  of  Flow 
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Fig.  9  Chamicol  Composition 
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Fig.  10  Vibrational  Energy 


38 


,  moles /gm  of  fluid 


AEDC-TDR-«4-29 


Fig.  n  I  ni  -  ni,  eg  I 
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